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Abstract 

A  numerical  model  is  being  constructed  to  simulate  the  movement  of  water  and 
constituents  within  a  coupled  groundwater  and  surface  water  system.  The  model  is  designed  to 
run  on  all  DoD  HPC  platforms.  Local  mesh  refinement  and  coarsening  and  domain- 
decomposition  preconditioners  have  been  integrated  into  the  model.  The  model  has  been 
exercised  successfully  on  several  test  problems.  Presently,  the  model  is  being  applied  to  an 
instrumented  watershed. 


Scientific  Problem 

The  Department  of  Defense  (DoD)  is  required  to  assess  the  environmental  impact  of  its 
activities  at  both  present  and  formerly  used  facilities.  When  warranted,  the  DoD  must  enact 
remedial  measures  to  address  environmental  problems.  The  potential  costs  associated  with 
environmental  remediation  at  DoD  sites  is  staggering.  In  addition  to  the  cost  of  remediation,  the 
DoD  risks  reduced  or  prohibited  access  to  its  training  facilities  unless  environmental  concerns 
are  addressed  adequately.  For  these  reasons,  accurate  environmental  assessments  and  effective 
remedial  designs  are  essential. 

Thorough  environmental  assessment  requires  that  the  ecosystem  be  examined  in  a  more 
holistic  fashion  than  is  customary.  Traditionally,  each  part  of  a  hydrologic  system  (groundwater, 
wetland,  drainage  basin,  etc.)  has  been  modeled  individually,  treating  the  other  parts  of  the 
hydrologic  system  as  sources  or  sinks  that  are  assumed  to  be  constant  or  described  with  simple 
empirical  functions.  Often,  the  disparate  temporal  or  spatial  scales  among  these  systems  justify 
this  uncoupled  treatment.  However,  in  some  cases,  the  degree  of  interaction  or  the  uncertainty 
in  the  magnitude  of  sources  or  sinks  requires  that  multiple  components  of  the  hydrologic  system 
be  considered  simultaneously.  Examples  include  groundwater-driven  hydrology  in  wedands, 
contaminant  exchange  between  surface  and  subsurface  systems,  and  heterogeneous,  transient 
infiltration.  In  many  locations,  such  as  south  Florida,  the  groundwater  and  surface  water 
systems  are  so  tightly  coupled  that  they  are  virtually  inseparable.  Figure  1  is  a  schematic 
showing  several  typical  points  of  interaction  between  groundwater  and  surface  water  systems. 
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Technical  Approach 


The  ADH  (ADaptive  Hydrology)  groundwater  model,  constructed  at  the  Engineering 
Research  Development  Center  (ERDC),  is  being  extended  for  use  as  the  model  foundation. 
Knowledge  gained  through  the  creation  of  the  HIVEL2D  surface-water  model,  also  at  the 
ERDC,  is  helping  develop  a  new  surface  water  modeling  capability  within  the  ADH  framework. 
The  new  model  is  named  SWGW. 

SWGW  couples  three-dimensional,  unsaturated  groundwater  modeling  to  two- 
dimensional,  shallow  water  modeling  at  the  surface.  Previous  attempts  to  account  for  the 
interaction  of  groundwater  and  surface  water  have  linked  a  saturated  groundwater  flow  model 
with  an  overland  flow  model  through  a  one-dimensional  (vertical)  representation  of  the 
unsaturated  zone.  Such  a  one-dimensional  unsaturated  flow  interface  typically  uses  a  simplified 
depiction  of  the  unsaturated  flow  as  wet-dry  interface  propagation  or,  even  simpler,  as  a 
homogenized-volume  represented  by  a  constitutive  equation.  The  SWGW  model  is  a  technical 
advance  because  it  maintains  full  three-dimensionality  in  the  unsaturated  zone,  permitting  the 
simulation  of  perched  aquifers,  inter-formation  flow  (lateral  flow  in  the  unsaturated  zone),  and 
infiltration  processes  in  heterogeneous  systems.  The  drawback  in  taking  this  approach  is  the 
large  computational  effort  required. 

Brief  Model  Description 

SWGW  approximates  the  solution  to  the  Richards  equation  for  groundwater  flow  and 
the  diffusive  wave  equation  for  surface  water  flow.  The  Richards  equation  is  a  combined  water 
balance  and  momentum  equation  for  saturated  and  partially  saturated  soil.  This  equation  is  non¬ 
linear  because  some  of  the  coefficients  (saturation  and  relative  permeability)  depend  on  the 
unknown  head.  The  diffusive  wave  equation  for  overland  flow  arises  by  neglecting  the 
acceleration  terms  in  the  full  St.  Venant  equations  (for  example,  Singh,  1996). 

Finite  elements  are  used  to  discretize  the  domain.  The  approximation  is  piecewise  linear 
in  space  and  piecewise  constant  in  time.  Groundwater  flow  is  solved  in  three  dimensions  using 
tetrahedra.  The  diffusive  wave  equation  is  approximated  on  triangles  that  comprise  a  surface  of 
the  three-dimensional  groundwater  flow  mesh.  Nodes  located  on  the  overland  flow  face  are 
dual  valued,  with  an  overland  flow  head  and  a  groundwater  head.  The  two  flow  regimes 
communicate  through  boundary  fluxes  computed  at  the  surface  of  the  groundwater  system.  By 
communicating  only  through  fluxes,  the  model  avoids  a  problem  found  in  other  models  coupling 
groundwater  and  surface  water.  Many  of  these  models  switch  boundary  conditions,  using 
Dirichlet  (head)  boundaries  when  the  depth  in  the  surface  water  is  non-zero,  and  Neumann 
(flux)  boundaries  for  recharge  when  the  surface  water  has  zero  depth.  Presently,  the  two  flow 
regimes  update  their  fluxes  only  at  each  time  step.  As  model  testing  proceeds,  it  may  be 
necessary  to  enforce  this  communication  at  each  non-linear  iteration. 

Mesh  Refinement/Coarsening 

Many  of  the  physical  problems  to  be  addressed  with  the  SWGW  model  contain  steep 
and  moving  spatial  gradients  in  the  solution  variables.  Examples  of  these  gradients  include  a 
moving  saturation  front  or  intermittent  well  in  the  groundwater  system,  a  traveling  wave  in  the 
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surface  water  system,  or  a  contamination  front  in  either  system,  for  example.  Capturing  these 
phenomena  with  a  fixed-mesh  model  would  require  extremely  fine  mesh  resolution  throughout 
the  domain.  Such  resolution  is  not  practical  for  many  problems  and  is  not  efficient  use  of 
resources  for  most  problems.  For  these  reasons,  the  SWGW  model  uses  local  mesh 
refinement  and  coarsening  to  add  and  delete  resolution,  as  needed,  to  capture  steep  gradients. 

Splitting  or  merging  elements  is  based  on  an  explicit  enor  indicator.  Presently,  the 
model  uses  an  inexpensive,  gradient-based  indicator,  but  more  accurate  (and  costly)  indicators 
(Schmidt,  1997)  are  available  in  the  model.  Elements  slated  for  refinement  are  divided  using  the 
edge  bisection  scheme  by  Liu  and  Joe  (1995)  (Figure  2).  Edges  are  ranked  by  the  refinement 
level  of  its  nodes  and  by  their  length.  The  ‘oldest’  edge  in  an  element  is  split  first.  If  all  elements 
are  the  same  ‘age’,  the  longest  edge  is  split  first. 

Parallelization 

SWGW  has  been  constructed  to  take  advantage  of  parallel  computer  architectures.  The 
domain  is  subdivided  in  a  data  parallel  scheme  shown  by  a  simple  example  in  Figure  3.  Nodes 
are  distributed  to  processors  uniquely.  Elements  that  fall  on  processor  boundaries  are  shared. 
Ghost  nodes  are  created  for  those  nodes  that  reside  off  processor,  but  contribute  to  a  shared 
element.  Border  nodes  are  on-processor,  but  are  seen  as  ghost  nodes  by  another  processor. 
Thus,  border  nodes  must  communicate  information  with  other  processors. 

Processor  partitions  generally  will  contain  elements  that  are  spatially  adjacent  to  each 
other  because  this  tends  to  minimize  inter-processor  communication.  Therefore,  as  the  mesh  is 
refined  locally  near  a  large  gradient,  a  majority  of  the  additional  elements  and  nodes  can  be 
created  on  only  a  few  processors.  Thus,  local  mesh  adaption  creates  an  inherent  workload 
imbalance  among  the  processors.  Periodic  repartitioning  is  needed  to  maintain  the  load  balance 
for  the  dynamic  system.  Mesh  refinement  occurring  on  one  processor  must  be  communicated 
to  other  affected  processors. 

Inter-processor  communication  is  handled  with  the  standard  Message  Passing  Interface 
(MPI)  libraries  to  ensure  the  model=s  portability  among  many  HPC  machines.  Thus  far,  the 
model  has  been  ran  on  the  Cray  T3E,  IBM  SP,  and  SGI  Origin  2000  at  the  Major  Shared 
Resource  Center  in  Vicksburg. 


Matrix  Preconditioners 

Recent  research  indicates  that  for  many  problems,  including  groundwater  transport  in 
naturally  heterogeneous  soils,  significant  resolution  is  necessary  to  produce  qualitatively  conect 
answers  (Tompson  and  Gelhar,  1990,  Howington  et  al,  1997).  This  revelation  comes  as  the 
trend  in  physical  problem  dimension  and  complexity  continues  to  increase  rapidly.  In  tandem, 
the  trend  toward  larger  physical  dimension  and  finer  resolution  is  leading  to  enormous  increases 
in  the  number  of  nodes  and  elements  in  a  typical  simulation.  The  SWGW  model  is  implicit  in 
time,  requiring  the  simultaneous  solution  of  large  non-linear  algebraic  systems.  The  number  of 
nodes  and  the  degrees  of  freedom  per  node  determined  the  size  of  this  system  of  equations.  An 
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inexact  Newton’s  method  is  used  to  linearize  the  problem.  Thus,  a  linear  system  must  be  solved 
for  each  Newton  iteration  and  several  Newton  iterations  may  be  required  for  each  time  step. 


A  general  concern  exists  with  the  solution  of  these  large  linear  systems.  The  number  of 
iterations  required  to  solve  the  system  grows  with  problem  size  for  most  schemes.  Therefore, 
the  work  required  for  matrix  assembly,  etc.  increases  linearly  with  the  number  of  processors, 
but  the  work  required  to  solve  the  system  can  grow  more  rapidly.  With  clever  preconditioning 
of  the  linear  system,  this  growth  in  the  number  of  iterations  can  be  dramatically  reduced 
(Tompson  et  al,  1994). 

A  domain-decomposition  approach  was  chosen  to  precondition  the  linear  system 
because  these  are  well  suited  for  parallel  implementation.  When  these  subdomains  are 
overlapping,  these  are  known  as  Schwarz  preconditioners.  The  preconditioner  options  in 
SWGW  are: 

•  Point  Jacobi 

•  One-level  Additive  Schwarz 

•  Two-level  Additive  Schwarz 

•  Two-level  Hybrid  Schwarz 

Point  Jacobi  preconditioning  makes  no  use  of  domain  decomposition.  The  remaining 
preconditioners  divide  the  domain  into  overlapping  subdomains.  Figure  4  shows  a  sample  fine 
mesh  and  four  subdomains. 

One-level  additive  Schwarz  is  simply  a  block  Jacobi  preconditioner.  A  fine-mesh  solve 
on  each  subdomain  is  followed  by  an  interpolation  back  to  the  full  preconditioning  matrix.  To 
extract  each  subdomain  from  the  larger  system,  one  must  assume  boundary  condtions.  Zero 
Dirichlet  boundary  conditions  are  used  on  the  subdomain  boundaries.  Two-level  additive 
Schwarz  schemes  add  a  full-domain  coarse  mesh  solve  to  the  subdomain  solves.  Basis  functions 
for  the  fine  mesh  elements  are  summed  to  create  a  single  basis  function  for  that  subdomain 
(Jenkins  et  al,  in  preparation).  This  summed  basis  function  for  one  of  the  four  subdomains  in  the 
sample  mesh  is  shown  in  Figure  5.  The  coarse  problem  consisting  of  a  single  matrix  entry  per 
subdomain  is  solved.  Parallelizing  the  coarse  mesh  is  not  yet  required  because  each  processor 
can  perform  the  coarse  solve  independently.  The  two-level  additive  and  two-level  hybrid 
schemes  combine  the  fine  and  coarse  mesh  information  differently.  This  domain  decomposition 
approach  is,  effectively,  a  simple,  multigrid  preconditioner  on  an  unstructured  mesh,  without  the 
complexity  of  creating,  maintaining,  and  parallelizing  multiple,  nested,  meshes. 


Application 

The  SWGW  model  has  been  applied  to  several  example  problems  and  application  to  a 
field  site  is  underway.  Among  the  example  problems  are  drainage  through  a  heterogeneous  soil 
column  and  rainfall/mnoff  in  a  simple  test  basin.  The  column  problem  is  intended  to  demonstrate 
the  model’s  capabilities  in  simulating  drainage  in  heterogeneous  soil.  The  simple  test  basin  is 
being  used  to  explore  the  fluxes  across  the  ground  surface  and  evaluate  the  performance  of  the 
overland  flow  model. 

A  field  problem  has  been  constructed  to  test  the  model.  The  Poplar  Creek  drainage 
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basin  at  Camp  Shelby  near  Hattiesburg,  MS  has  been  studied  extensively.  The  numerical  mesh 
to  be  used  in  the  initial  simulations  for  this  field  site  is  shown  in  Figure  6.  Rain  will  be  applied  to 
the  surface  and  flux  in  the  creek  will  be  compared  to  measured  discharge. 


Additional  Benefits 

This  work  was  motivated  by  the  need  to  simulate  the  interaction  between  groundwater 
and  surface  water.  However,  by  adopting  a  strategy  that  attaches  problem- specific  routines 
(groundwater  flow  equations  and  overland  flow  equations)  to  a  single,  computational  engine, 
several  software  development  and  maintenance  advantages  have  become  apparent.  The 
computational  engine  contains  matrix  solvers,  mesh  adaption  routines,  generic  finite  element 
routines,  and  parallel  communication  routines.  These  difficult  and  time-consuming  parts  of  the 
code  remain  virtually  unaffected  by  many  changes  to  the  groundwater  or  surface  water  routines. 


Another  major  issue  with  a  development  of  this  magnitude  is  code  maintenance.  If 
written  in  a  modular  fashion,  a  single  code  may  perform  simulations  for  groundwater,  surface 
water,  or  coupled  systems  with  little  overhead  penalty.  Therefore,  this  code  may  circumvent  the 
need  to  maintain  individual  codes,  which  often  contain  similar  modules.  By  sharing  related 
modules,  advancing  capabilities  are  kept  in  step  for  each  of  the  potential  applications.  Our  hope 
for  the  future  is  to  extend  the  model  to  include  other  components  of  the  hydrologic  system.  The 
difficulty  lies  in  keeping  the  code  at  a  manageable  size  and  in  keeping  the  components 
sufficiently  modular  to  permit  enhancements  and  maintenance. 


Status  and  Plans 

The  SWGW  model  presently  solves  coupled  groundwater  and  surface  water  flow  using 
a  diffusive  wave  approximation  for  overland  flow.  The  model  is  being  tested  against  simple 
problems  for  which  analytical  approximations  are  possible  and  against  data  from  instrumented 
watersheds.  Plans  exist  to  upgrade  the  surface  water  model  to  solve  the  full  shallow  water 
equations.  A  method  to  handle  canals  efficiently  in  one  dimension  will  be  explored.  A 
constituent  flux  will  be  added  at  the  boundary  between  the  flow  regimes  to  accommodate 
constituent  transport  between  surface  and  subsurface  waters. 

Because  the  mesh  partitioning  among  subdomains  (not  simply  among  processors) 
defines  the  coarse  mesh  for  the  Schwarz  preconditioner,  the  partitioning  scheme  deserves 
careful  study.  Preliminary  simulations  indicate  that  the  convergence  rate  is  sensitive  to  the  shape 
of  these  subdomains.  Likewise,  there  is  a  delicate  balance  when  choosing  the  number  of 
subdomains.  Having  too  many  subdomains  creates  a  very  large  coarse  mesh  problem,  while 
too  few  subdomains  requires  the  solution  of  large  problems  for  each  subdomain. 
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Figure  1.  Typical  interaction  of  groundwater  and  surface  water. 
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Figure  2.  Schematic  showing  the  refinement  and  coarsening  of  a  single  tetrahedral  element, 
adding  or  removing  one  element  and  one  node. 


Figure  3.  Schematic  showing  the  mesh  partitioning  scheme.  On  the  left,  a  simple  mesh  is 
divided  among  three  processors.  Nodes  are  assigned  to  processors.  Elements  may  be  shared 
by  processors.  On  the  right  is  the  view  from  processor  0  showing  interior  nodes,  border  nodes, 
and  ghost  nodes. 
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Figure  4.  An  example  fine  mesh  with  element  boundaries  shown  in  black  and  four 
preconditioner  subdomains.  The  light  gray  elements  are  the  overlapping  regions. 
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Figure  5.  The  coarse-mesh  basis  function  is  produced  by  summing  all  the  fine-mesh  basis 
functions  within  the  subdomain.  The  resulting  coarse-mesh  basis  function  is  constant  except  in 
the  elements  shared  with  the  other  subdomains. 
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Figure  6.  The  initial  Poplar  Creek  watershed  computational  mesh. 
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